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Abstract 

The symmetric heavy-light ansatz is a method for finding the ground state of any dilute un- 
polarized system of attractive two-component fermions. Operationally it can be viewed as a 
generalization of the Kohn-Sham equations in density functional theory applied to iV-body density 
correlations. While the original Hamiltonian has an exact Z2 symmetry, the heavy-light ansatz 
breaks this symmetry by skewing the mass ratio of the two components. In the limit where one 
component is infinitely heavy, the many-body problem can be solved in terms of single-particle 
orbitals. The original Z2 symmetry is recovered by enforcing Z2 symmetry as a constraint on 
iV-body density correlations for the two components. For the ID, 2D, and 3D attractive Hubbard 
models the method is in very good agreement with exact Lanczos calculations for few-body systems 
at arbitrary coupling. For the 3D attractive Hubbard model there is very good agreement with 
lattice Monte Carlo results for many-body systems in the limit of infinite scattering length. 



I. INTRODUCTION 



Dilute two-component fermions with attractive interactions have some universal features 
of relevance to several areas of physics. In three dimensions the unpolarized ground state 
is believed to be superfluid with properties in between a Bardeen-Cooper-Schrieffer (BCS) 
fermionic superfluid at weak coupling and a Bose-Einstein condensate (BEC) of bound 

n u n 

dimers at strong coupling 12J, |3J. The crossover transition occurs somewhere near the 
unitarity point, a scale- invariant point where the range of the interaction is zero and the 
scattering length is infinite. Much recent interest on this topic has been stimulated by 
experimental progress with cold atomic Fermi gases. Starting with a dilute Fermi gas, where 
the effective range of the interaction is negligible compared with the interparticle spacing, 



the scattering length can be tuned using a magnetic Feshbach resonance [4|, 

m 

In nuclear physics the phenomenology of the unitarity point is relevant to cold dilute neutron 
matter. The neutron scattering length is roughly —18 fm while the range of the interaction is 
comparable to the Compton wavelength of the pion, m^ 1 wl.4 fm. Therefore the unitarity 
point is approximately realized when the interparticle spacing is about 5 fm. This range of 
neutron density is expected in the inner crust of neutron stars. 

In this paper we introduce a new theoretical method called the symmetric heavy-light 
ansatz. The ansatz is used to find the ground state energy of dilute unpolarized two- 
component fermions with attractive interactions. While our main interest is the three- 
dimensional system, we demonstrate that the ansatz is accurate for any number of spatial 
dimensions. Operationally it can be described as an iV-body generalization of the Kohn- 
Sham equations in density functional theory. But rather than solving for one-body densities, 
we use the symmetric heavy-light ansatz to determine iV-body density correlations. While 
useful as a numerical technique, the ansatz also provides a heuristic picture of the underlying 
competition between Fermi repulsion and attractive interactions. 

We start with an effective Hamiltonian describing two-component fermions with short- 
range attractive interactions. For the unpolarized case the ground state has an exact Z2 
symmetry corresponding with interchanging components. The first step of the ansatz is to 
break this Z2 symmetry by skewing the mass ratio of the two components. In fact we take 
the extreme limit where one component is infinitely heavy. In this limit the many-body 
problem is completely solved in terms of single-particle orbitals. The original Z2 symmetry 
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is reintroduced as a constraint on the iV-body correlations for the two components. In the 
simple approximation where only the lowest orbitals are filled, this constraint completely 
determines the iV-body density correlations for each fermion component. We run the 
method through several numerical tests and find very good agreement for a number of exact 
Lanczos and Monte Carlo results for the attractive Hubbard model in one, two, and three 
dimensions. The low computational scaling of the ansatz allows us to make predictions for 
much larger systems which are not yet accessible using alternate methods. 

The organization of the paper is as follows. The beginning is a short summary of sym- 
metries of the effective Hamiltonian. After this we introduce the heavy-light Hamiltonian 
i^HL and the symmetric heavy-light ansatz. Some results relating to heavy-light orbitals are 
derived, and the spectrum of the heavy-light Hamiltonian is explored numerically for 
the unpolarized four-body system in three dimensions. This analysis provides motivation 
for the lowest filling approximation and the construction of iV-body fixed-point densities. 
We then design a Markov chain algorithm to generate fixed-point densities. As a first test 
we compare the results of the symmetric heavy-light ansatz with exact Lanczos results for 
the four-body system in three dimensions. We also compare with exact results for four- 
and six-body systems in both one and two dimensions. We then consider larger systems in 
three dimensions at unitarity and the ground state energy for general values of the scattering 
length. We conclude with a discussion of some extensions of the method and future work. 



II. SYMMETRIES OF THE HAMILTONIAN 



In continuum notation we can write the effective Hamiltonian for two-component fermions 
in d < 3 dimensions with short-range interactions as 




Our main interest is the case d — 3, but we also consider lower dimensional systems d = 1,2. 
CL/j and (V a are annihilation and creation operators for fermions with two components. We 
refer to these components as up and down spins. The mass of the fermion is m, and the 
coefficient C is assumed to be negative so that the interaction is attractive. The strength of 
C depends on the scheme used to regulate the short distance behavior. This Hamiltonian 
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has a global U(l) fermion-number symmetry, 



"a T (f)~ 




"a T (r)" 









(2) 



where is any real constant. It also has a global SU(2) spin symmetry 



(3) 

where a denotes the 2x2 Pauli spin matrices and (ft is any constant real three-component 
vector. Since there is no coupling between spin and orbital angular momentum, this SU(2) 
symmetry should be regarded as an internal symmetry decoupled from spatial rotations. 

The lowest-dimensional local bosonic operator that can be constructed from the annihi- 
lation field operators is 

ift 2 (r) = ai(r)ai(r). (4) 

We note that ift 2 is invariant under the SU(2) spin symmetry but phase rotates under the 
U(l) fermion-number symmetry, 



(5) 



Therefore if there is some critical temperature below which ift 2 has long-range spatial corre- 
lations, 

lim (^t( r > 2 (0)W0, (6) 

\r\— »oo \ / 

then the U(l) fermion-number symmetry is spontaneously broken. This condition of off- 



diagonal long-range order [111, [12| is the standard definition for superfluidity with S-wave 
pairing. We observe that the SU(2) spin symmetry is not broken by the expectation value 
of ip 2 . 



III. SYMMETRIC HEAVY-LIGHT ANSATZ 



Let and be the kinetic energy operators associated with the up and down spins 
respectively, 

1 



2m 
1 



dVaj(f)V 2 a T (r), 



-— / d d r a\(r)V\(r). 



(7) 
(8) 
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We define TThl as 

H mj = H-K,+K l = 2K l + V, (9) 

and refer to 77hl as the heavy-light Hamiltonian. In 77hl we have deleted the kinetic energy 
for the up spin while doubling the kinetic energy of the down spin. This could be viewed 
as introducing spin-dependent masses m-\ and mj_. For the original Hamiltonian H we have 
= mi = m, while in 77hl we have m\ = oo, mi = m/2. 

The physics of the two-body system in the center-of-mass frame is exactly the same for 
H and TThl- The reduced mass \x defined by 

i = -U-L do, 

[i mi mi 

equals m/2 in both cases. The exact equivalence of H and Hbl for the center-of-mass two- 
body system is preserved by most regularization schemes such as dimensional regularization, 
momentum cutoff schemes, and Hamiltonian lattice regularization. Hamiltonian lattice 
regularization is discussed in the appendix. 

While i?HL and H are identical for the two-body system, they are very different for more 
than two particles. For example when the ultraviolet cutoff scale goes to infinity for the 
three-dimensional system, H has a well-defined continuum limit while Hul is unbounded 
below due to a clustering instability. The simplest example where this occurs is the three- 
body system consisting of two up spins and one down spin. It is known that this system 
collapses due to an attractive 1/r 2 potential when > 13 mi [3]. On the other hand the 
clustering instability in 77hl is eliminated when we project onto quantum states invariant 
under the interchange of up and down spins. In the following we show that the minimum 
expectation value for 77 H l restricted to this up-down symmetric space equals the ground 
state energy of H. 

We denote a state with iVf up spins and iVj_ down spins as an N^,Ni state. We also 
specify the total momentum P and total spin S of the SU(2) spin representation. Let 
I ^ at at) be the normalized ground state of H for the N, N system in a periodic cube of 
length L. Since the 577(2) spin symmetry is not broken by the expectation value of i/j 2 , we 
assume that I^jvtv) lies in the spin-invariant sector 5 = 0. Let E% N be the corresponding 
ground state energy. Let U be any unitary operator which interchanges the up and down 
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spins through a 7r-radian spin rotation. Without loss of generality we take U to be 



U = exp 

Clearly 



%- \ d d r a\(r)a^(r) - I d d f a\{f) ai {f) 



(11) 



rfc^U = -ia u (12) 

U^a t U = -ia T . (13) 

The phases appearing in Eq. ffl2l and ffl3|) are necessary in order that the pair annihilation 
operator ifj 2 (r) remains invariant, 

U^ 2 (r)U = rfa^a^U = -a^a^r) = ^ 2 (f). (14) 

Any state with an even number of fermions is invariant under two successive transforma- 
tions of U. Therefore U acting on the space of even fermion states generates a unitary 
representation of the group Z 2 . 

Let |$jv,7v) be any N, N state which is Z 2 invariant, 

U\<£> N , N ) = |$jv,tv) • (15) 

Since 

(&n,n\ \<& N) n) = (<5>n,n\ K l \® N>N ) , (16) 

it follows that 

($jv,iv| H nL \& N<N ) = (<5>n,n\ H |$jv )JV ) . (17) 

In particular this means 

(®N,N | -^HL | *&N,n) ■ {$N,n\ H \&n,n) /10 x 

mm — — ; = mm — — — ; — . 18 



Since is a ^-invariant state we conclude that 

(®N,N #HL \®N,n) r-iO , ln % 

mir ! , — ik ^ \ — = E n,n- (19) 



U\&N,N / = \ V N,N j \-*-JV,JV | ■■■ iV,iV / 

We refer to this exact relation as the symmetric heavy-light ansatz. The task of computing 
the ground state energy E° N N is reduced to finding the minimum of value of the Rayleigh-Ritz 
ratio 

(&n,n\ Hul\&n,n) ^2q^ 
under the constraint of Z 2 invariance, 

U\$ N , N ) = \$n,n)- (21) 
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IV. HEAVY-LIGHT ORBITAL S 



In this section we discuss some properties of the energy eigenstates of the heavy-light 
Hamiltonian -£/hl- Since the up spins are infinitely massive, we can fix these particles at 
locations Ri, R 2 , ■ • • , Rn- Due to antisymmetry each Ri must be distinct, and we use the 
shorthand notation 



R 



(22) 



for the unordered set of vectors. We define |R) as the corresponding antisymmetric state 
of localized particles, 

|R) = -^=YjSgn{ir) R^ <8> R^ 2) ) ® ■ ■ ■ R\{n)) ■ (23) 

' 7T 

The summation is over all permutations tt of the integers 1, 2, ■ • • ,N. sgn(ir) equals —1 for 
odd permutations and +1 for even permutations. The normalization for |R) is then simply 

(R'|R) =detA ii (R',R), (24) 

where 

Ay (R', R) = 5^ (jl'i — R^j . (25) 

Later in the discussion when we consider lattice models the Dirac delta functions will be 
replaced by Kronecker delta functions. 

For any given R the down spins see a static delta function potential from each up spin, 

C8^(f-Ri). (26) 

i=l,«- ,N 

For the very special case where R is a regular cubic array, this system is known as the 
Kronig- Penney model [ijj]. In general though R is irregular without translational and 
rotational symmetries. For a given R let the normalized single-particle orbitals be |/ 3 -(R)) 
with corresponding eigenvalues Ej(H). By convention we label the orbitals so that Ej(H) 
increases with j, 

Ei(R) < E 2 (R) < < Ej(R) < ■ ■ ■. (27) 
For a single orbital we write the position-space wavefunction as 

/,(f,R) = (f|/,(R)). (28) 
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Analogous to R, we define 



r = {ri,f 2 ,--- ,fjv} 



(29) 



for the unordered set of distinct vectors fj. These correspond with possible locations for 
the N down spins. We also define the antisymmetric product of position eigenstates 



with normalization 



-L= E s 9n(n) |rv ( i)) <g> |fv (2) ) <g> • • • \f n{N) ) . 



(r' |r) = detA^r). 



(30) 



(31) 



Since there is no interaction between down spins, each eigenstate of Hwl for fixed R is 
an antisymmetric product of single-particle orbitals. Let us write the normalized quantum 
state with N down spins filling orbitals ji, • • • ,jw as 

f\ \f jn (R)) = -4j E s ^^) lAoi)( R )> ® |Ao 2 )(R-)> ® • • • . (32) 



n=l,-,N 

The position state wavefunction for this state is a Slater determinant, 



(r| A l/ J n(R))=det 



/^(f^R) / j2 (fi,R) ••• / jjv (fi,R) 
/^R) R) ••• / JJV (f 2 ,R) 

fji(r N ,R) f j2 (f N ,R) ■■■ f jN (f u R) 



(33) 



All eigenstates of Hkl are a tensor product of up-spin position eigenstates and down-spin 
orbitals, 

n=l,-,N 

with corresponding eigenvalue 

E ^( R )- ( 35 ) 

n=l,- ,7V 

We can use the eigenstates of Hhl to perform a basis decomposition of the ground state 
I^jvjv)- For any given R and orbital indices j±, ■ ■ ■ , ]n, we write the inner product with 
the ground state as 



<R|® A (/?'«( R ) 

n=l,-,JV 



(36) 
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with normalization 



We find that 



[d dN Rj2\F* 

J {Jn} 



{in} | 



(37) 



E 



N,N 



(*° N>N \ H\*° NtN ) = (K,n\ Hul \K,n) 
f d dN R^\Fn, {jn} \ 2 £ E Jn (R). 

J {in} n=l,- ,iV 



(38) 



m I I ^ 

The coefficients LFr,^} define a normalized probability distribution. Weighted by this 
probability distribution, E° N N is the average over all R and {j n } of the orbital energy sums 

E Ej n (R). (39) 

n=l,- ,iV 

The energy constraint Eq. (|38p is an exact relation satisfied by the ground state energy 
E° N N of the original Hamiltonian H. Since the phase of -Fa,{j n } is irrelevant, solving this 
constraint is qualitatively easier than solving for the ground state wavefunction l^jviv)- I n 
this respect Eq. (1381) can be regarded as the starting point for a generalized density functional 
approach. However instead of single-particle densities we work with many-body densities 
given by the coefficients -Fr, {.,•„} . 



V. H AND H nL ON THE LATTICE 



Throughout this section we discuss Hamiltonian lattice regularization for interacting two- 



component fermions with short-range interactions, 
and Euclidean lattice formulations can be found in 
26. 
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urther details 
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27] . 



13, 



19 



or 



Doth Hami 
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21 



22 



23 



tonian 
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A. Attractive Hubbard model 



On a (i-dimensional spatial lattice we can write H as 



H = Kt + K\ + V, 



where 



^- ^ [2a{ (r)a t (r) - a\{f)a^{f+ I) - a\(r)a^(r - t) 

r,l 



(40) 



(41) 



r,l 

and 



^I = ^E [ 2a K*0«i(*0 - a\(f) ai (r + 1) - a\(f) ai (f- /)] , (42) 



V = cY J ^M^ a ^ a ^- ( 43 ) 



Here f is an integer-valued <i-dimensional spatial lattice vector, and I = l,...,d are lattice 
unit vectors in each of the spatial directions. We write to for the fermion mass and C for 
the coupling constant. Throughout we use dimensionless parameters and operators, which 
correspond with physical values multiplied by the appropriate power of the spatial lattice 
spacing a. This lattice model is the same as the ci-dimensional attractive Hubbard model. 
The Hamiltonian for the Hubbard model is usually written as 

^Hubbard = ~t Yl [4(rK(r + 1) + 4 (fK(f - f)] + U a\{r)a\(r)a^r) ai (r). (44) 
Therefore in d dimensions, 

^Hubbard =H + {N, + - (45) 

and 

U = C, t = ±. (46) 
We define the heavy-light Hamiltonian on the lattice as 

H HL = H-Ki + K l = 2K l + V. (47) 

As in the continuum case the change from if to if hl can be viewed as altering the masses of 
the spins. In if we have toj = mj = ra, while in if hl we have raj = oo, toj_ = to/2. More 
generally we can define 

if (0) = ff + (-fsT T + K l )9 = (l- 0)tf T + (1 + 0)fq + V, (48) 

which interpolates between H at = and ff hl at = 1. In terms of the spin-dependent 

masses toj and raj., 

= ™T- m * . (49) 

TO| + TO j 
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FIG. 1: Low-energy spectrum of H{6) at total momentum P = for the four-body system on 
a 4 3 lattice The energies are written as a fraction of the ground state energy E®'^ for the 
non- interacting system. The coupling constant C is tuned to the unitarity point. 



B. Spectrum of the four-body system in three dimensions 

To gain some intuition for the low-energy spectrum of H(6), we consider the unpolarized 
four-body system, iV-j- = N± = 2, in three dimensions. We tune the coupling constant C 
to the unitarity limit where the scattering length is infinite. This procedure uses Liischer's 
finite- volume scattering formula and is discussed in the appendix. In Fig. [T] we show the 
low-energy spectrum of H(9) for the four-body system at total momentum P = on a 



4 lattice. The energies are computed using the Lanczos method [28] and written as a 



fraction of the ground state energy E^^ for the non-interacting system. The concentric 
circles indicate exact degeneracies as well as approximate degeneracies too close to visually 
distinguish in the plot. 

The ground state energy remains relatively constant for 9 < 0.6. For larger 9 the ground 
state energy decreases more substantially while the density of low-energy states increases. 
The data at 9 = 1 corresponds with the spectrum of the heavy-light Hamiltonian i^HL- From 
Fig. [1] we count 35 energy states at 9 = 1 with E 2 ^/ E^^ less than 0.25. These lowest 
35 energy states correspond with the 35 independent states obtained by placing down spins 
in the lowest two energy orbitals for each possible up-spin configuration R = {R\, R2}, and 
then projecting onto P — 0. A sketch of the one-body density distribution of down spins 
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FIG. 2: Sketch of the one-body density distribution of down spins in the lowest orbitals for a fixed 
configuration of up spins. 

in the lowest orbitals for a fixed configuration of up spins is shown in Fig. [2j 

VI. LOWEST FILLING APPROXIMATION 

A. Lowest orbital filling and effective distribution |-Fr| 2 

On the lattice we write a discrete sum over R rather than a continuous integral. Therefore 
Eq. (|38|) becomes 

e n,n = J2J2\ f *>im\ 2 E ^«( r )- ( 5 °) 

R {j„} n=l,-,N 

In our example of the four-body system in three dimensions at unitarity, we found that for 
each R the sum 

71=1,2 

is numerically close to E® 2 only when the down spins occupy the lowest orbitals, {j n } = 
{1,2}. Given that the weighted average 

EEI F ^ } | 2 E^( R ) ( 52 ) 

R {j n } n=l,2 
2 

equals E% 2 , we conclude that -Fk,{j n } is dominated by the set of lowest orbitals {j n } = 
{1,2} for each R. This suggests a general approximation scheme which keeps only the 
lowest orbitals. 
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As in our four-body example at unitarity, let us assume that most of the weight of the 

i 1 2 

normalized distribution ■Pk,{j„} lies in orbital sets {j n } where all orbitals lie near the 
bottom N in energy. Let 

F * = jEl^if- ( 53 ) 

V 

Then 

E °n,n = J2J2\ F ^}\ 2 E E jn (R)^J2\ F ^\ 2 E (54) 

R {j n } n=l,-,N R j=l,-, N 

This approximation uses an effective distribution with weight \Fn\ which is nonzero only 
for {j n } = {!,•••, N}. We refer to this as the lowest filling approximation. 



B. iV-body fixed-point densities 

In the lowest filling approximation we are restricted to states of the general form 



R 



|R>® A l/i( R )> 

3=1,-, N 



(55) 



Ei^i 2 = L ( 56 ) 



R 



In general any |$at,at) of this form will not satisfy the .^-invariance condition U \<&n,n) = 
\&n,n) exactly. Therefore we search for a weaker symmetry constraint that follows from Z 2 
invariance but can be satisfied exactly in the lowest filling approximation. In the following 
we design a constraint which is exactly solvable in the lowest filling approximation, with a 
unique solution |-F R | 2 for each R. 

Let us define the up-spin and down-spin particle densities, 

p T (f) = a\(r)a r (r), (57) 

Pl(r) = (58) 
Let G[pi,pi\ be any functional involving the particle densities. Then any state |3>jv,jv) 
satisfying the invariance condition U \&n,n) = \®n,n) also satisfies 

($jv,Ar| : G [p h Pi] : \^n,n) = (®n,n\ ■ G [p h p T ] : \<5>n,n) . (59) 
In particular we have 

($n,n\ : Pi(n) x ■ ■ ■ x pi(r N ) : \$ n ,n) = ($n,n\ ■ P;(fi) x • • • x pi(r N ) : \$ n ,n) • (60) 
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The : symbols denote normal ordering. Upon summation this constraint implies analogous 
relations for the j-body densities for each j < N, 

(®n,n\ ■ Pt( f i) x • • • x p T (r}) : \$ N , N ) = (<$>n,n\ ■ Pi(n) x • • • x p^r}) : \<S> N , N ) . (61) 
We can write the iV-body up-spin densities in terms of F R , 



(®n,n\ ■ Pi(Ri) x • • • x Pj (-Rat) : \$n,n) 

= ($ NlN \ al(R N ) x • • • x aUR^a^Rx) x • • • x a^(R N ) \®n,n) = N\ |F R | 



(62) 



We can write the iV-body down-spin densities in terms of F R and the inner product 
<r| A l/;(R)>, 

3=1,-, N 



($n,n\ ■ Pi(n) x • • • x Pl (f N ) : \$ N , N ) = N\Y^ 



R 



j=h-, N 



x |Fi 



R 



(63) 



If the iV-body densities for up spins and down spins are equal then 



R 



(r| A ^'( R )> 

3=1,-, N 



x|F R | 2 = |F r | 2 . 



(64) 



We refer to any |F R | which satisfies this relation an iV-body fixed-point density. The 
existence and uniqueness of fixed-point solutions |F R | will be proved in the next section. 

At this point we comment on the clustering instability of i^HL- In two dimensions when 
any two of the up-spin locations Ri and Rj come close together, the lowest heavy-light orbital 
energy in physical units has a divergence proportional to a -1 , where a is the lattice spacing. 
In three dimensions the divergence is a" 2 . However any iV-body fixed-point distribution 
|F R | must vanish as 
in 

(r| A l^( R )>- 

3=1,-, N 

This is enough to remove the clustering instability for d < 3. Therefore the continuum limit 
is well defined for the lowest filling approximation to the ground state energy, 



Ri — Rj 



due to antisymmetry with respect to exchange of r-i and fj 

(65) 



J? 



3=1,-, N 



(66) 



R 



2 

where |F R | is an iV-body fixed-point density. 
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C. Lowest filling approximation for the four-body system in three dimensions 



We test how well the lowest filling approximation works for the unpolarized four-body 
system in three dimensions. Again the coupling is tuned to the unitarity point and the 
lattice volume is 4 3 . Once we project onto the total momentum P = subspace, there 
are 35 independent configurations R for the up spins. Similarly there are 35 independent 
configurations for the down spin positions r. We construct a 35 x 35 matrix M (r, R) with 
elements 



M(r,R) 



(67) 



(r| A IA'( R )> 

J"=l,2 

where |/i(R)) , |/ 2 (R)) are the lowest orbitals for R. We then solve for the iV-body fixed- 



point density |FrJ satisfying 



^M(r,R) |F R | 2 = |F r | 2 . 



(68) 



R 



With |F R | the ground state energy can be calculated using 

R j=l,2 



(69) 



We can also compute two-body correlation functions. The same-spin correlation function 
in the lowest filling approximation is 

<Wf) = £X*2,2| : p^f+f^p^f') : |*° 2 > 

r' 

= ^<^ i2 |:p T (f+f')p T (r-"):|^ 2 > 



E { S Ri,i 
R 



s 



(70) 



In our dimensionless lattice normalization the same-spin correlation function summed over 
r equals 2. Also the opposite-spin correlation function is given by 



r' 

= E<*kl +?>,(?'): |*° 2 ) 



EE 

R .7=1,2 



/,(f+ J R 1 ,R) + /,(r + J R 2 ,R) 



IF, 



R 



(71) 
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FIG. 3: Comparison of exact and heavy-light results for the same-spin and opposite-spin correlation 
functions. The data is binned by radial distance in lattice units rounded to the nearest half integer. 

The opposite-spin correlation function summed over f equals 4. Fig. [3] shows a comparison 
of exact results and heavy-light results in the lowest filling approximation for the same-spin 
and opposite-spin correlation functions. The exact results are Lanczos calculations for the 
ground state of H. The data is binned by radial distance r = |r| in lattice units rounded 
to the nearest half integer. The agreement between exact and heavy-light results for the 
correlation functions is rather good. The agreement for the ground state energy is also 
quite good. The exact value for the ratio E^/E®'^ is 0.138 while the heavy-light result 
gives 0.128. This level of quantitative agreement may seem surprising given the simplicity 
of the lowest filling approximation to the heavy-light ansatz. The method simply patches 
together a collection of single-particle orbitals using the iV-body fixed-point constraint in 
Eq. ( 1641) . However the method appears to accurately describe the competition between 
the short-range attractive interaction and Fermi repulsion. We now check whether the 
agreement is also good for different interaction strengths. 

The ratio E^ n/E^ 1 ^ tends towards negative infinity at strong attractive coupling and 
is therefore somewhat unwieldy to plot over a large range of coupling strengths. Let us 
define a different dimensionless ratio which is convenient for the attractive Hubbard model 
at arbitrary coupling and arbitrary dimensions. At fixed lattice volume L d we define the 
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ratio 



T-,0 77,0, free 



— iv, im — iv,iv l<-?c\\ 

6N ' N ~ N IE \ + E ^ ( } 

iV 1-^1,1 1 ^ -^jv.jv 

where E^ is the ground state energy for the two-body dimer in the same lattice volume L d . 
In the limit of weak attractive coupling, e^,N tends towards from below. In this case the 
physics is dominated by Fermi repulsion, and the ground state is similar to that of a free 
Fermi gas. In the limit of strong attractive coupling, cn,n tends towards —1 from above. 
Here the ground state consists of tightly-bound dimers with only weak interactions between 
dimers. For e^jy in the midrange between and —1, the competition between attractive 
interaction and Fermi repulsion is stalemated to some intermediate balance point. The 
unitarity limit in three dimensions fits this category. In the limit of large iV at unitarity 
the N xj term can be neglected relative to E^ r ^ e . In this case 

/?0 77iO,free 

^N,N N,N _ F Q /p 0,frce , / 7 oN 

N,N 

In Fig. H] we show results for e2,2 for lattice lengths L = 3, 4, 5 and various couplings. We 
express the coupling as a ratio of Hubbard model parameters 

U/t = 2mC. (74) 

The plot on the left shows exact results computed using the Lanczos method. The unitarity 
point corresponds with U/t = —7.914, and the signal of scale invariance at unitarity can be 
seen by the agreement in e2,2 for different L. The plot on the right shows the difference 
between heavy-light results in the lowest filling approximation and exact Lanczos results. 
From weak coupling at e2,2 ~ to strong coupling at e2,2 ~ — 1, the error of the heavy-light 
calculation is bounded by 0.02. The method appears to accurately describe the crossover 
from four weakly-interacting fermions to a condensed pair of bosonic dimers. 

VII. MARKOV CHAIN FOR iV-BODY FIXED-POINT DENSITIES 

In this section we design a Markov chain process which generates the iV-body fixed-point 
density |-Fr| 2 that solves the constraint in Eq. fl64l) . Our construction establishes both 
existence and uniqueness. Let us define 



T(r,R) 



(r| A l/i( R )> 

j=l,-,N 
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(75) 
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FIG. 4: Comparison of four-body exact Lanczos results, e| X 2 Ct , and heavy-light results, ei? \i f° r 
L = 3, 4, 5 lattices in three dimensions plotted versus coupling U/t = 2mC. 

We note that T(r, R) > for every pair r, R. Also for each R, 

£T(r,R) = l. 

r 

We define a Markov chain 

r(o) ^ R(D ^ > R (*) ^ R (*+U _> . . . (77) 

with transition probability 

p (R (fc+1) | R (fc) ) = T(R( fc+1 ), R«). (78) 

Each state in the Markov chain is a set of N distinct <i-dimensional vectors on an L d 
lattice. We note that 

T(R( H1 >,R W ) = (79) 
if and only if R^ fe+1 ^ and f\ /j(R")| are exactly orthogonal. Since the state space 
is finite the possibility of accidental orthogonality between |R^ +1 ^ and f\ |/j(Rw)^ 

J=1,-,JV 

is a set of measure zero and can avoided by an arbitrarily small change in the coupling 
C. However it is also possible that the orthogonality arises from some mismatch of exactly 
conserved quantum numbers. These quantum numbers would be associated with some sym- 
metry subgroup of the lattice shared by R^ fc+1 ) and RW, It might be a reflection symmetry, 
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(76) 



rotational symmetry, translational symmetry, or some combination of each. However for 
large L the relative proportion of such symmetric configurations R/ fc+1 ) and is sup- 
pressed by powers of L and therefore exceedingly rare. We let e be a small positive number 
and define a modified chain with transition probability matrix 

p(R* +1 > R<*>) =T [ (R< t+1 »,RW), (80) 
T W = "^, . (81) 

r 

For any R' and R" there exists some chain of finite length from R' to R" such that the 
transition probability at each step is nonzero. In fact we can get there in only one step. 
Hence the Markov chain is ergodic and there exists a unique invariant distribution, call it 
|F R | 2 , such that 

^T e (R,R)x|F R | 2 = |F R ,| 2 . (82) 

R 

2 

We take the limit as e — > and obtain \Fr\ as the unique iV-body fixed-point density. 
This fixed-point Markov chain shares some features with iterative solutions of the Kohn- 



Sham equations in density functional theory 29J . Both involve finding the lowest orbitals of 



a Schrodinger potential which in turn depends on the particle densities. However there are 
important differences between this method and the Kohn-Sham equations. The fixed-point 
Markov chain presented above solves for iV-body densities rather than one-body densities. 
Also in contrast with density functional theory, this method yields an ab initio solution 
for the uniform system. We discuss possible extensions to non-uniform systems in the 
discussion section. 

As a precise test of the Markov chain algorithm for \Fr.\ we compare with heavy-light 
calculations presented earlier for the unpolarized four-body system at unitarity on a 4 3 
lattice. For this calculation and all other heavy-light calculations presented in the following 
we generate the fixed-point |Fr| using 4000 steps of the Markov chain 

R(0) _> R(D ^ „ R(fc) ^ R (fc+D _> . . . ( 83 ) 

with transition probability 

p (R {k+1) | R (fe) ) = T(R( fe+1 ), R(*>). (84) 
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FIG. 5: Comparison of Markov chain and direct heavy- light calculations for the same-spin and 
opposite-spin correlation functions 

Each step of the chain RS k ' — > R( fc+1 ) is produced using a Metropolis algorithm with 100 
updates per up-spin particle location according to the probability distribution 



TfR, R (fc) ^ 



(R| A lM R(fc) )> 

3=1,-, N 



2 



(85) 



The entire calculation is repeated several times with different random number seeds and the 
results are averaged. The standard deviation of the distribution is used to determine the 
stochastic error of the average. 

In Fig. [5] we show a comparison of calculations for the same-spin and opposite-spin cor- 
relations functions using the Markov chain algorithm and the direct calculations shown 
previously in Fig. [3j The error bars shown are estimated stochastic errors. We see that 
the Markov chain algorithm reproduces the direct calculations for the correlation functions 
with no detectable systematic error. The Markov chain heavy-light result for E\ 2 jE, 



0,free 
2,2 



is 0.127(1), which agrees with the direct heavy-light calculation of 0.128. The advantage 
of the Markov chain algorithm is that the computational scaling for large systems is very 
favorable, nearly linear in the lattice volume and number of particles. 
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VIII. FEW-BODY SYSTEMS IN ONE AND TWO DIMENSIONS 



In this section we compare heavy-light calculations for ejv,jv m the lowest filling approxi- 
mation with Lanczos calculations for few-body systems in one and two spatial dimensions. 
The rapid computational scaling of the Lanczos calculation, L^ -1 ^, limits the system sizes 
for which the comparison is possible. Let us first discuss the properties of bound dimers in 
the ID, 2D, and 3D attractive Hubbard models. 

In one dimension at infinite volume the dimer has binding energy {^J 

B 2 = -El = ^ + ..., (86) 

and characteristic size 

1 2 

7^ = ^\c\ + " ' ■ (87) 

The ellipses denote lattice spacing corrections which become important when the character- 
istic size is not much larger than one lattice spacing. In two dimensions at infinite volume 
the dimer has binding energy [3^] 

^ n 3.24tt 2 / 4tt \ 

and characteristic size 

1 1 ( 2tt \ 

exp + • ■ ■ . (89) 



yJmB~ 2 1.80tt V m C 

n t 



3, 



iree dimensions at infinite volume and mC < —3.957, the dimer has binding energy 
231 
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B E° 1 167r2 ( 1 1 V 

1,1 wascatt m \mC 3.957/ 1 " 

and characteristic size 

1 1 

= M^ + A) + "" (91) 

In one and two dimensions the attractive Hubbard model has no nontrivial scale-invariant 
point analogous to the unitarity point in three dimensions. 

Fig. shows e2,2 for the unpolarized four-body system in one dimension. Fig. [7] shows 
e3 t 3 for the unpolarized six-body system in one dimension. For both the four-body and 
six-body systems the error of the heavy-light result is smaller than 0.015. We note that the 
six-body heavy-light data for L = 4 is exact. This is because the system can be regarded 
as the ground state of one up-spin hole and one down-spin hole. Similar to the unpolarized 
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FIG. 6: Comparison of four-body exact Lanczos results, e| X 2 Ct , and heavy-light results, e^\, for 
L = 4, 6, 8, 10, 12 lattices in one dimension plotted versus coupling U/t = 2mC. 
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FIG. 7: Comparison of six-body exact Lanczos results, e™ *, and heavy-light results, e^, for 
L = 4, 6, 8, 10, 12 lattices in one dimension plotted versus coupling U/t = 2mC. 



two-particle system, the ground state of the unpolarized two-hole system depends only on 
the reduced mass /i, 

1 1 1 

(92) 



1 1 1 



which is unaltered in the heavy-light formalism. 

Just as in the three-dimensional example at unitarity, the method appears to accurately 
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FIG. 8: Comparison of four-body exact Lanczos results, e| X 2 Ct , and heavy-light results, e^JjS for 
L = 3,4, 5, 6, 8, 10 lattices in two dimensions plotted versus coupling U/t = 2mC. 

describe the crossover from fermions to condensed bosonic dimers in one dimension. The 
difference between the ground state energy E Q N N and N times the dimer energy E^i can 
be viewed as Fermi repulsion of overlapping of dimer wavefunctions. In one dimension the 
dimer wavefunction has an exponential tail proportional to 

exp i-raC \x\ \ = exp ( — \x\ J . (93) 



2 1 7 1 \4t 

This simple picture is consistent with the exponential dependence of cn,n on U/t seen in 
both the four-body and six-body results for large negative U/t. 

Fig. [8] shows e2,2 for the unpolarized four-body system in two dimensions. Fig. [9] shows 
e3 ; 3 for the unpolarized six-body system in two dimensions. In both cases the error of the 
heavy-light result is bounded by 0.02. Just as in the one- and three-dimensional examples, 
the method appears to accurately describe the crossover from fermions to condensed bosonic 
dimers. In two dimensions the dimer wavefunction at infinite volume has an exponential 
tail with characteristic length 

1 1 / 2tt \ 1 / 4tt \ 

ex P = 7^T ex P "7777 • ( 94 ) 



1.80tt rnC J 1.80vr U/t t 

The dimer increases its size dramatically for small negative U/t. But at finite volume it 
eventually wraps around the periodic boundary at distance L. This boundary effect changes 
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FIG. 9: Comparison of six-body exact Lanczos results, e™ 01 , and heavy-light results, eljg, for 
L = 3,4 lattices in two dimensions plotted versus coupling U/t = 2mC. 

the behavior of e^^N- This can be seen as a crossover to approximately linear dependence 
on U/t for U/t~> —4 in both the four-body and six-body results. 

IX. DIMER-DIMER SCATTERING IN THREE DIMENSIONS 

We have already presented results for e2,2 in three dimensions. In this section we re- 
analyze the data to extract low-energy scattering parameters for dimer-dimer scattering. 
Luscher's formula 
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32 



34| relates the energy levels for any two-body system in a 
finite periodic cube to the S-wave phase shift. In the appendix we discuss how to use 
Luscher's formula in the two-body system with one up spin and one down spin to measure 
the fermion-fermion scattering length, a sca tt- At strong coupling we can also use Luscher's 
formula in the four-body system to measure the S-wave dimer-dimer phase shift. When L 
is much bigger than a scatt we can interpret the energy difference 

AE° 2>2 = E° 2j2 - 2E° ltl (95) 

as the energy of the two-dimer system relative to threshold. We then determine the dimer- 
dimer phase shift using 

>v 



p D cot5 DD = -^rS(rj) , V=\-^-) , (96) 



24 



a 
a 

CO 

o 



Q 

a 

ft, 



0.5 



-0.5 - 
-1 - 
-1.5 



-2.5 





1 1 1 

Exact Lanczos 


A 




Heavy-light ansatz 
a DD - 0.60 a scatt 








A A _ 


A 
A 


A 

A A A A 




- 








1 1 1 





0.5 



(PD°scatt) 



1 1.5 

2 



FIG. 10: Results for pd cot <5dd as a function of the dimer momentum squared. Both quantities are 
measured in units of a sca tt • For comparison we also show the limit of pb cot <5dd at zero momentum 
determined by the result odd ~ 0.60a scat t. 



where is the dimer momentum, S (rj) is the three-dimensional zeta function defined in 
the appendix, and <5dd is the dimer-dimer S-wave phase shift. It is convenient to measure 
everything in units of the fermion-fermion scattering length a sca tt- 

Results for prj cot <5dd are shown in Fig.fTUlas a function of dimer momentum squared. For 



comparison we also indicate the result odd ~ 0.60a scatt [35j, [36( which determines pv cotfcD 
at zero momentum using the effective range expansion, 

p D cot <5 DD w — + ~tddPd H • (97) 

odd * 

We see in Fig. [TU] that the heavy-light and exact Lanczos results agree at the level of a 
few percent level. Also both agree with the dimer-dimer scattering length result odd ~ 
0.60a scatt . From this plot we can also estimate the dimer-dimer effective range, 

r DD w 2.6a scatt . (98) 
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FIG. 11: Heavy-light results for £nn for N from 2 to 32. The results are plotted versus L , 
where L ranges from 4 to 16. 

X. BCS-BEC CROSSOVER IN THREE DIMENSIONS 



A. Many-body results at unitarity 

We present symmetric heavy-light results in the lowest filling approximation at unitarity 



for 



6v, 



N 



^ N,N 
r-i0,free 
N.N 



(99) 



for a wide range of values for N and L. Fig. [TT] shows results for £n,n for the four-body 
system, iV = 2, up to the sixty- four-body system, N = 32. The results are plotted versus 
where L ranges from 4 to 16. The dependence on L is relatively mild. The maximum 
for £,n,n for N = 7 appears to be caused by the closed shell at N = 7 for the free Fermi 
system in a periodic cube. In the continuum and thermodynamic limits L, N — > oo we find 



£= lim ^ = 0.31(1). 

L,N— *oo 



(100) 



The heavy-light results can be compared with published Euclidean lattice Monte Carlo 
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results at small volumes shown in Table 1 
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Table 1: £at,jv using Euclidean lattice Monte Carlo 



N = 3 


N = 5 


N=7 


N = 9 


N =11 


0.19(2) 


0.24(2) 


0.28(2) 


0.23(2) 


0.25(2) 



These results correspond with an average of data from L = 4, 5, 6 for each TV = 3, 5, 7, 9 
and L = 5,6 for iV = 11. For the same N and L, the heavy-light results in Fig. [TT1 agree 
reasonably well with the Euclidean Monte Carlo data, with a relative difference of at most 
20%. 

The heavy-light data for larger values of L probe much further towards the continuum 
limit. These results suggest that lattice cutoff effects cannot explain the discrepancy between 
lattice Monte Carlo results and continuum fixed-node Monte Carlo results for the same values 
of N. Fixed-node Green's function Monte Carlo calculations have found £jv,iv to be 0.44(1) 



37J,|38( and 0.42(1) [39(] for comparable values of N. Further studies are needed to determine 



if different nodal surfaces can produce lower ground state energies in fixed-node calculations. 
The symmetric heavy-light ansatz may be useful in probing this question more deeply. In 

there are a 

number of other theoretical calculations 
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531 ] and experimental measurements of £ 



40JJ41L [43, 143| 

i y, y, y, w 
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45 



46 



48 



49 



50 



51 



54j which span the range from 



about 0.2 to 0.6. 



B. Results for general scattering length 

We consider £jv,jv as a function of kpa scaitt . In the limit of strong attractive coupling and 
the thermodynamic limit iV — ► oo at fixed density, we get 

p0 iup0 , N(N-l) 47ra Dp 

where E\ x is the energy for one dimer and odd is the dimer-dimer scattering length. The 
first term takes into account the binding energy of the dimer while the second gives the 
contribution due to dimer-dimer interactions. Although we have written the expansion 
in powers of fcF^scatt, a more accurate estimate of the appropriate expansion parameter is 
°scatt divided by the average spacing between particles d ~ (67r 2 ) 1 / 3 /c^ 1 . We expect the 

°(^F a scatt) err ° r to be Sma11 for ^Fflscatt < 2. 
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For a periodic cube we have 

1 / 9 iV\ 2/3 (6tc 2 N) 2/3 , s 

e ' = *r[**v) < 102 > 

and 

*,= <*^. (103) 
In the continuum limit the energy for one dimer is 

Ki = \— . ( 104 ) 



ma 2 



scatt 



and the dimer-dimer scattering length is approximately odd ~ 0.60a scat t- Putting these 
together we get 

E° 5 5 

^ = l in L = " op „2 + - 60 X T^-^ascatt + 0(4<&att). ( 105 ) 

JV 00 0/t F a scatt i07r 

We refer to the first two terms in this expansion as leading order (LO) and next-to-leading 
order (NLO) in the strong-coupling expansion. 



In the weak-coupling limit we have 
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K,n 10 , 4(11 -2 In 2) 

+ 9tt ^ scatt + 21vr2 



= = 1 + ^T^catt + ^7T73p^ fe F a Ltt + °(4 a Ltt)- ( 106 ) 

<J N,N 

We refer to the first three terms in this expansion as leading order (LO), next-to- leading 
order (NLO), and next-to-next-to- leading order (NNLO) in the weak-coupling expansion. 
In Fig. [12] we show as a function of fc^ 1 a~^ tt for iV = 32 and L = 16 using the symmetric 
heavy-light ansatz in the lowest filling approximation. For comparison we show the analytic 
strong-coupling and weak-coupling results. 

We see that the heavy-light results are very close to the strong-coupling results for 
^p la scatt ^ 0-3- The lowest filling approximation is also not bad in the weak-coupling 
limit for kp a~ catt ~ — !■ Near the unitarity point we can expand 

E° N 

€ N > N = p 0,free = £ _ ^ lk F °Wtt ~ °^att ^ ■ ( 107 ) 

N,N 

For iV = 32 and L = 16 we find by least squares fitting that 

£ = 0.306(1), (108) 
5 = 0.805(2), (109) 
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FIG. 12: Plot of £n,n f° r N = 32 and L = 16 as a function of fcp & SC att- ^or comparison we show 
the analytic strong-coupling and weak-coupling results. 



6 = 0.63(3). (110) 
The heavy-light result for £i is within 20% of a Euclidean lattice Monte Carlo calculation 



25j | which found £i = 1.0(1). In contrast with £, there is genera 



literature on the value of calculated using different methods [38 



agreement in the recent 
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XI. DISCUSSION 



A. Lowest filling approximation and beyond 

In this analysis we have tested the symmetric heavy-light ansatz in the lowest filling 
approximation on numerous few- and many-body systems for the ID, 2D, and 3D attractive 
Hubbard models. In each case the method appears to be accurate. It is perhaps not 
surprising that the lowest filling approximation is reliable at strong attractive coupling. In 
this limit the fermions form tightly-bound dimers which are weakly interacting and condense 
in the ground state. In the ground state the relative momentum between dimer pairs is 
relatively low, and the up-spin coordinates R can be viewed as surrogates for the locations of 
dimer pairs. The iV-body fixed-point constraint sets |.Fr,| so that up spins and down spins 
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are arranged in the same manner relative to neighboring particles. At very strong coupling 
|Fr| is essentially independent of R, indicating that the dimers are non-interacting bosons 
each in the zero momentum state. As the attractive coupling becomes weaker, the size of 
the dimers increases and the effect of Fermi repulsion becomes more important. 

For the systems we have considered here the lowest filling approximation works well at 
weak attractive coupling. However in future studies one may consider systematic improve- 
ments to the lowest filling approximation that include excited orbitals. This improvement 
would probably be needed to study the repulsive Hubbard model where there is no gap 
between the lowest N orbitals and higher orbitals. In the lowest filling approximation the 
coefficients I-FrJ 2 are set by imposing Z 2 symmetry on the same-spin correlations, 

($n,n\ ■ Pt(n) x • • • x Pl(r N ) : \& N , N ) = ($n,n\ ■ Pi(n) x • • • x Pl (f N ) : \& N , N ) . (Ill) 

If excited orbitals are included, a larger set of coefficients Fr^j must be determined. 
These can set by imposing Z 2 symmetry on mixed-spin correlations such as 

{&n,n\ ■ Pt(^i) x Pt(^) x Pi(f 3 ) : \$n,n) 

= ($ N , N \ : Pi(fi) x Pl (r 2 ) x p T (f 3 ) : \& N , N ) . (112) 

It is not clear that an algorithm can be designed to solve these constraints a priori, as we 
did for the lowest filling approximation. Therefore one approach would be to determine 
some parameterization for Fr,^} with unknown coefficients setting the relative weight of 
each orbital set {j n } for fixed R. These coefficients could then be determined a posteriori 
by least squares fitting to the mixed-spin correlation constraints. 

B. R-commuting operators and phases 

We define an R-commuting operator as any operator which commutes with the up-spin 
density P]{r) for all r. Some examples of R-commuting operators include iJjjL and the 
up-spin and down-spin density operators. For any R-commuting operator O the ground 
state expectation value is 

(Kn\0\K,n) = J2 E F l{u} F ^}Of Jnhm , (H3) 
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where 



U {Jn},{j' n } 



(R\®/\(f)M) 



o 



|R)®AI^( R )> 



(114) 



Suppose O is an R-commuting M-body operator with M <^ N. Due to the orthogonality 
of orbitals, the orbital sets {j n } and {j' n } must be the same for all orbitals left untouched 
by O. Because of this constraint the diagonal elements O^-^ ^ are enhanced by powers of 
N relative to the off-diagonal matrix elements r jf y. For the special case O = Hrl all 

the off-diagonal matrix elements are in fact zero. If the weight of -pR,{j„} | is dominated by 
orbital sets {j n } close to the lowest orbital filling {1, • • • , N}, we can approximate O^-^ 
by the diagonal element at lowest orbital filling, 



U {jn},{jn} ~ U n,- 



{1,- ,N}, {1, -,N}- 



(115) 



This leaves us with 



{in} 



,AT},{1,- ,JV} l-^R 



(116) 



R 



R 



The expectation values of R-commuting operators are simple since we can keep only 
diagonal matrix elements, and the sign and phase of is irrelevant. The expectation 
values of other operators are more challenging. One example that is not R-commuting is 
the difermion pair correlation, 

^ 2 W(0), (117) 



where 



^ 2 (f) = a^rja^r). 



(118) 



Here we need to compute matrix elements for orbitals from different up-spin configurations 
as well as nontrivial geometric phases which may arise going from one up-spin configuration 
R to another R'. 



C. Non-uniform systems 

The Markov chain algorithm for iV-body fixed-point densities at lowest filling works also 
for the non-uniform case with external potential V(r). Heavy-light calculations should be 
possible for systems such as harmonic traps used in cold atomic experiments. In the case of 
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slowly-varying potentials these calculations could be compared with results obtained using 



58 



59 



60 



61 



62|. 



density functional theory 

Let us define H [V] with potential V(f) as 

H [V] = H + Vtf Mf) + Pl (r)} , (119) 

f 

where H is the attractive Hubbard Hamiltonian defined in Eq. (HOI) . If N [V] is the ground 
state energy of H [V], the straightforward generalization of the symmetric heavy-light ansatz 
gives 

min { * N *l H ™M\*™) =] % N M, (120) 

u\$ N , N ) = \<f> NtN ) \™N,N \™N,N) 

where 

Hul [V] = H HL + J2 V(f) [p T (?) + Pl(?] ■ (121) 
f 

However we could also define a more general form for i^HL, 

H HL [V, V A ] = H HL + J2 [V(? - V A (f)} p T (f) + W(?) + Va(?) P|(r), (122) 

r r 

for arbitrary Vj^if). The contribution of Va(?) cancels from the expectation value in 
Eq. (I120p . For unpolarized systems adding an overall constant to the auxiliary potential 
Va(t) has no effect, and so we can set 

E v *(?) = °- ( 123 ) 

The optimal V A {? can be found by minimizing the Rayleigh-Ritz ratio, 



mm ; — : ; — , (124) 

C/|*jv,jv) = |*jv,jv) \&N,N \®N,N/ 

in the lowest filling approximation. Roughly speaking the Va(?) adjusts the single-particle 
down-spin density, 

($ N)N \ : Pl {r) : \$ N>N ) , (125) 

while the iV-body fixed-point constraint insures that the single-particle densities are equal 
for both spins, 

,n\ '■ P^{? '■ \&n,n) — ($n,n\ '■ Pi?) '■ \&n,n) ■ (126) 

Applications of this approach to non-uniform systems will be discussed in a future publica- 
tion. 
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XII. SUMMARY 



We have presented a many-body approach called the symmetric heavy-light ansatz. It 
is an approximate method for finding the ground state energy of dilute unpolarized two- 
component fermions with attractive interactions. Although the Hamiltonian has an exact 
Z2 symmetry, the ansatz breaks the symmetry by changing the ratio of the masses of the 
two components. We considered the extreme limit where one component is infinitely heavy 
and the many-body problem can be solved in terms of single-particle orbitals. The original 
Z 2 symmetry was reintroduced by setting the A^-body density correlations for the two com- 
ponents equal for all N. A Markov chain algorithm was designed to generate exactly this 
A^-body density constraint. 

To test the method we first compared the results of the symmetric heavy-light ansatz 
with exact Lanczos results for the four-body system in three dimensions. We then tested 
the method with exact results for four- and six-body systems in both one and two dimen- 
sions. We then considered dimer-dimer scattering and larger systems in three dimensions 
at unitarity and arbitrary values for the scattering length. The results indicate that the 
method is quite accurate and robust. We discussed some extensions beyond the lowest 
filling approximation and applications to non-uniform systems. 



XIII. ACKNOWLEDGEMENTS 



The author thanks Gautam Rupak and Thomas Schafer for discussions and comments. 
This work is supported in part by DOE grant DE-FG02-03ER41260. 



APPENDIX A: FINITE- VOLUME SCATTERING IN THREE DIMENSIONS 



In this appendix we use Luscher's formula 
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32 



33 



34[ to relate the coefficient C in the 



three-dimensional Hubbard model to the S-wave scattering length. We consider one up-spin 
particle and one down-spin particle in a periodic cube of length L. Luscher's formula relates 
the two-particle energy levels in the center-of-mass frame to the S-wave phase shift, 

2 



j9cot<5 (p) = —S (77) 

7lL 



Lp 



(Al) 
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where S(r]) is the three-dimensional zeta function, 

9(A 2 - n 2 ) 



S{rj) = lim 

A^oo 



E 



47rA 



n z — rj 



For \rj\ < 1 we can expand in powers of rj, 



1 



S(r)) = h lim 



1 



E 



fl(A 2 - ft 2 ) 
n 2 — rj 



AttA 



- + S + S^ 1 + S 2 v 2 + S 3V 3 ■ ■ ■ 
V 



where 



So = lim 



A^oo 



E 



d(A 2 - n 2 ) 



- 47rA 



n- 



Si = V — ^TT j > 1. 



The first few coefficients are 



For small momenta the effective range expansion gives 

f f 



p cot 5 (p) 



+ 7T r 0P + 



C^scatt 2 

where a sca tt is the scattering length and tq is the effective range. 
In terms of rj, the energy of the two-body scattering state is 



P 2 1] /27T 

m m\L 



(A2) 



(A3) 

(A4) 
(A5) 



S = -8.913631, Sx = 16.532288, S 2 = 8.401924, S 3 = 6.945808, 

S 4 = 6.426119, S 5 = 6.202149, S 6 = 6.098184, S 7 = 6.048263. (A6) 



(A7) 



(A8) 



We compute the two-particle scattering pole at energy -E po i e by summing the two-particle 
bubble diagrams shown in Fig. [13j For spin-dependent masses m-f and mj_ the single-particle 
dispersion relations are 



(p) = — ( X ~~ cospi "> ' 
wi(p) = — V(l-cosp z ). 



(A9) 
(A10) 
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; = l ;>o-o<; 

n = 0,l,... 1 n 



FIG. 13: Sum of bubble diagrams contributing to two-particle scattering 

After summing the bubble diagrams as a geometric series, the relation between C and E vo \ e 
is 

- i = lim -L V i ^ . (All) 

C L ^°° L 5 ^ + ^(2vrfc/L) + u^nk/L) 

When combined with Eq. (1Alj) . (1A7|) . and fjA8j) . we can relate C and the scattering length 

°scatt- We note that for both H and H HL we have 

2 

^(2TxkjL) + u^nk/L) = — [1 - cos (27r^/L)] . (A12) 

m ; 

Therefore the coefficient C is exactly the same in both cases. The unitarity limit corresponds 
with the value U/t = 2mC = -7.914. 
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